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This  paper  addresses  the  problem  of  multicriteria  (versus  single-criterion)  parametrical  identification  of 
the  autonomously  controlled  cargo  parafoil.  Based  on  the  structural  identification  as  an  initial  step  toward 
creation  of  an  adequate  model  of  the  parafoil,  a  high-fidelity  model  including  several  dozens  of  optimization 
parameters  has  been  developed.  The  present  paper  proposes  the  correct  statement  of  the  multicriteria 
parametrical  identification  problem  including  the  necessity  to  investigate  the  feasible  set  of  variable 
parameters.  The  paper  advocates  the  use  of  the  Parameter  Space  Investigation  method  and  Multicriteria 
Optimization  /  Vector  Identification  software  package  to  solve  the  problem. 


I.  Introduction 

AS  well  known,  the  problem  of  the  development  of  the  mathematical  model  of  some  dynamic  object  includes 
two  necessary  stages.  First,  the  equations  of  motions  governing  system’s  dynamics  should  be  derived  from  the 
nature  of  the  system.  Once  this  first  step  in  constructing  a  mathematical  model,  the  so-called  structural  identification 
(defining  a  number  and  type  of  equations  of  motions),  has  been  completed,  the  next  step,  the  known  as  a 
parametrical  identification,  i.e.  finding  numerical  values  of  variable  parameters  to  better  match  experimental  data 
should  be  carried  on. 

While  the  structural  identification  for  parachute-  and  parafoil-based  payload  delivery  systems  is  considered  to 
be  more  or  less  settled,1"5  the  parametrical  identification  (defining  aerodynamic  and  control  coefficients,  apparent- 
mass-tensor  elements,  etc.),  especially  for  those  high-degree-of-freedom  models  developed  in  the  past  decade1,3,5 
still  needs  to  be  addressed  and  it  is  being  addressed  by  different  group  of  researches  for  different  aerodynamic 
deceleration  systems. 

By  their  nature,  applied  identification  problems  are  multicriteria  problems.  However,  as  a  rule  these  problems 
have  been  treated  as  single-criterion  problems.6,7  Usually  it  is  done  by  using  the  most  important  criterion,  or  by 
using  several  criteria,  but  one  at  a  time.  The  standard  approach  however  is  to  develop  a  single  compound  criterion 
that  weighs  criteria  relative  to  their  importance. 

By  present  time,  several  single-criterion  approaches  have  been  developed  and  used  to  identify  the  parameters  of 
different  payload  delivery  systems.  Kurashova  and  Vishnyak8  used  maximum  likelihood  method  to  determine  the 
aerodynamic  characteristics  of  gliding  parachute.  They  suggested  identification  of  longitudinal  parameters  using 
experimental  data  obtained  from  lorry  equipped  with  attachment  points,  measurement  and  control  system.  Their 
method  exhibited  verification  errors  of  less  then  10-15  percent. 

Jann9  discussed  application  of  system  identification  methods  (maximum  likelihood  as  well)  to  the  acquired 
database  of  the  parafoil-load  system  (ALEX-I).  Thereafter,  he  determined  the  essential  parameters  of  the 
autonomous  landing  system.  He  described  how  the  incorporated  parameters  were  estimated  and  discussed  the  results 
and  their  applicability.  He  developed  two  different  mathematical  models  which  describe  the  real  system.  One  was  a 
3 -DoF  model  and  another  -  a  4-DoF  model.  However  capabilities  of  these  were  limited  as  they  do  not  account  for 
the  distance  between  center  of  mass  and  aerodynamic  reference  point.  Later,  these  models  accounted  for  the 
actuators  which  move  the  control  lines.  Jann  also  presented  an  approach  for  theoretical  calculations  of  the 
aerodynamic  coefficients  based  on  the  extended  lifting  theory  and  validated  those  using  real  flight  test  data  on 
powered  parafoil  ALEX10. 

Kothandaraman  and  Rotea11  developed  a  SPSA  (Simultaneous  Perturbation  Stochastic  Approximation) 
algorithm  for  parameter  estimation  used  for  nonlinear  parachute  model.  The  SPSA  is  a  tool  for  optimization  that 
doesn’t  rely  on  a  costly  gradient  computation.  They  claimed  their  method  is  useful  where  many  parameters  are  to  be 
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optimized.  They  used  this  method  to  determine  three  aerodynamic  coefficients,  four  apparent  mass  coefficients,  and 
the  initial  states  of  the  G-12  parachute  model. 

Rogers12  used  an  extended  Kalman  filter  algorithm  implementation  to  estimate  the  aerodynamic,  wind,  mass 
property  and  measurement  errors  for  controlled  low-glide  parachutes.  This  implementation  incorporated  two  new 
approaches:  i)  attitude  error  formulation,  to  eliminate  the  mathematical  singularity  associated  with  vertical  flight, 
and  ii)  incorporation  of  apparent  masses  as  a  part  of  motion  dynamics.  Rodgers  presented  details  of  the  linearization 
of  the  nonlinear  equations  of  motion  and  measurement  equations,  as  well  as  a  summary  of  the  system  error 
dynamics.  His  results  based  on  simulated  data  showed  that  aerodynamic  characteristics  and  winds  can  be  estimated 
separately  from  apparent  mass  coefficients. 

Hur  and  Valasek13  also  chose  the  Observer/Kalman  Filter  Identification  (OKID)  methodology  for  identification 
of  the  longitudinal  and  lateral/directional  dynamical  models  of  the  Buckeye  parafoil- vehicle  system.  OKID  is  a  time 
domain  technique  which  identifies  a  discrete  input/output  mapping  from  known  input  output  data  records.  Since  first 
being  developed  by  Juang  in  the  early  1990’s,14  the  method  has  been  successfully  employed  to  identify  linear  system 
models  of  flexible  spacecraft  structures  and  aircraft.  The  dynamics  of  the  Buckeye  vehicle  were  modeled  with  8 
DoF:  six  for  the  parafoil,  and  two  for  the  relative  pitch  and  yaw  attitudes  of  the  vehicle.  Based  on  preliminary  results 
the  authors  drew  a  conclusion  that  the  OKID  method  can  identify  the  parafoil-vehicle  dynamic  system  effectively 
and  accurately. 

Analysis  of  these  attempts  that  address  model  verification  leads  to  the  conclusion  that  major  differences 
between  them  lie  in  the  way  the  authors  account  for  the  influence  of  numerous  interrelated  parameters  on  the  motion 
of  entire  system.  The  present  paper  does  not  intend  to  discuss  various  numerical  techniques  of  parameter 
identification,  but  addresses  the  physical  issues  (multicriteria  essence)  behind  the  identification  process  and  is 
organized  as  follows.  Section  2  introduces  the  cargo  parafoil  model  developed  during  a  structural  identification5  and 
names  several  groups  of  parameters  to  be  identified.  Because  of  the  different  nature  of  these  groups  of  parameters 
several  adequacy  criteria  are  suggested  to  be  used.  That  raises  the  necessity  to  employ  a  multicriteria  parametrical 
identification  technique.  First,  Sections  3  and  4  formulate  multicriteria  optimization  problem  and  introduce  the  PSI 
method  developed  to  manage  such  problems.  Then,  Section  5  shows  how  multicriteria  optimization  routine  can  be 
converted  to  a  multicriteria  identification  algorithm.  The  paper  continues  with  Section  6  where  corresponding 
software  is  introduced.  While  an  extensive  identification  experiment  is  still  continues,  the  results  of  preliminary  runs 
are  discussed  in  Section  7. 


II.  Structural  versus  Parametrical  Identification 


For  some  generic  parafoil-payload  system  Pegasus,  the  issue  of  parametrical  identification  has  already  been 
addressed  earlier.5  Two  additional  degrees  of  freedom  (payload  pitching  and  yawing)  added  to  the  existing  set  of  six 
differential  equations  yields  8-DoF  model  that  can  be  concisely  represented  in  the  following  form: 


’  F  ’ 

' 

◄ 

II 

E-i 

•V 

-LAV* 

> 

(1) 

1 

M,_ 

> 

K  =R£(ac)oc 

■  A, 

=  RE(Ap)(ap 

(2) 

;rv, 


(3) 


where  vector  V*  =  [vr,o^,co^ J  consists  of  a  velocity  vector  of  the  system’s  center  of  gravity  V  =  [u,v,wf  , 
angular  velocity  vectors  of  parafoil  (canopy)  (oc  =  \pc,qc,rcf  and  payload  wp  -  [o,  qp ,  rp  ]  ,  vectors  Ac  and  Ac 

(A  =  [(p,0,y/f )  constitute  Euler  angle  triads  for  the  canopy  and  payload  respectively,  P.  presents  an  inertial 

position  of  the  system.  On  the  right-hand  side  of  Eq.(l),  block  3-by-3  matrix  A  (where  each  element  is  a  3-by-3 
matrix  itself)  represents  a  mass-inertia  matrix  of  the  system,  vectors  F,  Mc  and  constitute  aerodynamic  force  and 
moment  acting  on  the  system  and  its  two  components  (canopy  and  payload  respectively),  and  matrix  £  is  a  block  3- 
by-3  matrix  where  each  element  is  a  3-by-3  matrix  too.  In  Eqs.(2),(3),  notations  R^(A)  and  stand  for  the 
matrix  operator  acting  on  a  vector  A  and  rotation  matrix  respectively.  The  reference  values  of  the  major 
aerodynamic  and  control  coefficients  are  determined  by  the  wind-tunnel  data  and  other  previous  studies.  Reference 
values  for  parameters  characterizing  the  bundle  and  some  other  uncertain  variables  are  also  added. 
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What  is  of  most  importance  from  the  standpoint  of  identification  is  that  some  of  these  vectors  and  matrices 
depend  parametrically  on  different  sets  of  variable  parameters.5  Yet,  since  the  coefficients  in  these  series  are  known 
with  some  probability  only,  it  is  suggested  to  multiply  these  reference  values  by  gain  factors.  Thus  tuning  the  model 
means  adjusting  values  of  these  gains  within  the  certain  (feasible)  range. 

By  analyzing  the  model,  sets  of  variable  parameters  were  gathered  into  the  following  groups:  i)  a  group  of 
aerodynamic  force  coefficients  kF  ;  ii)  a  group  of  moment  coefficients  kM ;  iii)  a  group  of  mass-geometry 

coefficients  kG ;  and  iv)  a  group  of  apparent  mass  coefficients  ka  .  For  instance,  the  aerodynamic-force-coefficients 
vector  kF  includes  following  coefficients:  ka<< ,  kCLa,  kCLS/,  kCDo,  kCDo,  kAj,  kASf,  kCY/jnom,  kCY/}grad,  kCYS^nom, 
kCY5  grad  j  group  of  mass-geometry  coefficients  kM  contains  krig  ,  kx<G ,  etc. 

Another  feature  is  that  during  real  drops  the  wind  profile  cannot  be  measured  simultaneously  with  the  rest  of  the 
states,  therefore  creating  an  obvious  uncertainty  for  the  identification  algorithm.  This  uncertainty  together  with  some 
other  uncertainties  in  the  system  geometry  and  state  of  control  surfaces  (not  all  state  variables  were  observed  and/or 
are  available)  produce  one  more  set  of  variable  parameters  k v  .  In  total  as  much  as  several  dozen  variable 
parameters  derived  from  the  model  (l)-(3)  need  to  be  identified. 

Consider  now  the  adequacy  criteria  that  can  be  used.  Of  course  they  are  completely  based  on  the  flight  test  data 
available.  To  this  end  the  flight  test  data  acquired  at  different  (from  4 Hz  to  100 Hz)  rate  by  the  global  positioning 
system  receiver  and  inertial  measurement  unit  installed  atop  of  payload  contains  the  following  information: 


•  local  tangent  plane  coordinates; 

•  components  of  inertial  velocity; 

•  attitude  of  the  measurement  unit; 

•  angular  rates; 

•  state  of  the  control  surfaces  (flaps). 

As  mentioned  the  uncontrolled  dropsonde  released  together  with  the  parafoil  gathers  current  wind  profile 
(horizontal  speed  components  versus  altitude). 

Analyzing  this  available  data  one  can  think  of  the  following  adequacy  criteria: 

•  proximity  of  simulated  trajectory  to  the  real  one  (that  can  be  further  split  into  the  horizontal  and  vertical 
components); 

•  closeness  of  the  speed/heading  profiles; 

•  adequacy  of  the  natural  eigenvalues  for  all  channels; 

•  closeness  of  system  response  to  control  actions. 

Therefore,  rather  than  using  a  single-criterion  identification  as  it  was  done  by  other  researchers  it  is  proposed  to 
use  a  multiple-criterion  identification  with  several  sets  of  different  parameters  grouped  by  their  influence  onto  the 
general  parafoil-based  cargo  system  performance  (see  Fig.  1). 


Input  ^  Inputs 

Figure  1.  Single-criteria  identification  (a)  and  multi-criteria  identification  (b). 


Necessity  of  employing  multicriteria  (or  vector)  identification  technique  can  be  understood  from  the  following: 

1.  Although  the  developed  model  of  a  cargo  parafoil  system5  seems  to  work  fairly  well  and  reflect  all  major 
features  of  the  real  object  dynamics,  generally  we  cannot  assert  a  sufficient  correspondence  between  the  model 
and  the  object.  This  obviously  limits  utility  of  the  single-criterion  identification  to  evaluate  adequacy  of  the 
model.  In  multicriteria  identification  there  is  no  necessity  of  artificially  introducing  a  single  criterion  to  the 
detriment  of  the  physical  essence  of  the  problem; 
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2.  The  richness  of  the  flight  test  data  allows  using  several  particular  proximity  criteria  to  evaluate  adequacy  of  the 
mathematical  model,  i.e.  determining  to  what  extent  the  mathematical  model  corresponds  to  the  physical 
system  in  principle; 

3.  The  fact  is  that  there  is  not  enough  preliminary  information  about  lower  and  upper  limits  for  many  of  the 
variables  meaning  that  prior  to  finding  the  optimal  ones,  we  should  be  interested  in  determining  the  feasible  set 
of  these  parameters,  as  well  as  performing  sensitivity  analysis. 

The  multicriteria  identification  is  a  relatively  new  direction  that  is  of  great  value  in  modem  engineering 
applications.15  The  numerical  technique  to  solve  such  kind  of  problems  has  usually  been  adopted  from  multicriteria 
optimization.  Multicriteria  optimization  methods  have  been  considered  in  many  articles,  monographs  and 
handbooks.16"18  However,  experts  continue  to  experience  difficulties  in  correctly  stating  optimization  problems  in 
engineering.  These  troubles  typically  emerge  when  trying  to  define  the  set  of  feasible  solutions,  i.e.  the  constraints 
imposed  on  the  design  variables,  functional  relationships,  and  criteria.  The  Parameter  Space  Investigation  (PSI) 
method19  was  developed  specifically  for  the  correct  statement  and  solution  of  engineering  optimization  problems. 
The  PSI  method  has  already  been  used  successfully  for  the  statement  and  solution  of  the  different  types  of 
multicriteria  problems  such  as  design,  design  with  control,  optional  development  of  prototypes,  finite  element 
models,  and  decomposition  and  aggregation  of  large-scale  systems.  It  was  also  implemented  for  identification  of  the 
static  systems.  Naturally,  we  would  like  to  employ  this  method  for  the  problem  at  hand. 

The  following  briefly  describes  the  essence  of  the  PSI-method  developed  initially  for  multicriteria  optimization 
problems  and  shows  how  it  can  be  used  for  identification  problem  at  hand. 

III.  Formulation  of  Multicriteria  Optimization  Problem 

Notice,  the  model  (l)-(3)  can  be  reduced  to: 

x  =  f(x,a,0,  te[t0;tf],  x,=(b  =x0.  (4) 

In  this  model  x  =  {xl9...9xn}  is  a  state  vector  with  initial  conditions  x0 ,  and  a  =  {al9...,a  }  is  a  vector  of  variable 
parameters  to  be  optimized. 

Consider  next  three  types  of  constraints  one  should  account  for  in  order  to  formulate  a  multicriteria  optimization 
problem  correctly.  They  are:  i)  parametric,  ii)  functional,  and  iii)  criteria  constraints. 

The  parametric  constraints  in  general  have  the  form 

a.<aj  <  aj ,  j  =  l,...,p.  (5) 

(For  mechanical  systems  a  .  usually  represent  geometrical  dimensions,  stiffnesses,  masses,  moments  of  inertia, 
damping  factors,  etc.,  and  define  a  parallelepiped  II  in  the  ^-dimensional  space.) 

The  functional  constraints  may  be  written  in  the  similar  form  as 

Cj  ^  fj(  a)  ^  Cj  ,  j  =  1  for  t  e  [t0;tf] .  (6) 

The  third  group  of  constraints  involves  local  quality  criteria 

/  =  l,...,v  (7) 

that  should  be  minimized/maximized.  Because  of  the  multicriteria  nature  of  the  problem  to  decrease  the  total 
number  of  the  reasonable  candidate  solutions  (avoid  the  situations  when  the  values  of  certain  criteria  are 
unacceptable  from  the  expert's  standpoint)  the  criterial  constraints  must  be  introduced 

$»<$;*,  /  =  l,...,v.  (8) 

Here  Op  is  the  worst  value  of  a  particular  criterion  expert  can  tolerate  while  ameliorating  other  criteria.  (Without 
loss  of  generality  here  and  further  on  we  consider  a  minimum  problem.) 

The  functional  dependences  /.  (a)  and  the  quality  criteria  Oz  (a)  may  be  functionals  of  the  interval  curves  of 
the  analyzed  differential  equations  (or  alternative  mathematical  models)  or  just  functions  of  a  . 

The  major  difference  between  criterial  and  hard  functional  constraints  is  that  the  values  of  Op  are  not  known 
beforehand  and  have  to  be  determined  while  solving  the  problem.  They  are  subject  to  expert's  revision  (he  either 
tightens  or  loosens  them).  For  the  sake  of  flexibility  the  functional  constraints  can  also  be  represented  in  the  form  of 
pseudocriteria,  especially  when  they  are  not  firm. 
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Obviously,  two  last  types  of  constraints  limit  initial  space  II  to  subspace  G,  G  c  IT  and  finally  to  some  feasible 
set  D,  DcGcII  as  is  shown  on  Fig.2. 

We  are  now  ready  to  formulate  the  multicriteria  optimization  problem  for  system  (4)  and  sets  of  constraints  (5), 
(6),  and  (8). 

The  multicriteria  optimization  problem  is  to  find  an  Edgeworth-Pareto  set  EPS,  EPS  cD,  so  that  the 
following  holds: 

<D/(EPS)  =  min<D/(a),  /  =  l,...,v  .  (9) 

aeD 

After  finding  EPS  the  most  preferable  or  optimal  vector  a0 ,  a0  e  EPS  can  be  finally  determined  (chosen). 

At  this  point  is  important  to  point  out  the  following.  Unlike  well-conditioned  traditional  single-criteria 
optimization  we  are  not  only  interesting  in  finding  a0  but  in  defining  the  feasible  and  Edgeworth-Pareto  sets  first. 


Figure  2.  Illustration  of  subsets  II,  G,  and  D  in  2D  space. 


IV.  Essence  of  the  PSI-Method 

The  PSI-method  is  based  on  populating  the  search  region  II  with  a  uniformly  distributed  sequence  of  points.  To 
produce  such  sequence  a  set  of  auxiliary  uniformly-distributed  on  the  unit  ^-dimensional  cube  points  Qu  i  =  1,...,M 
is  generated  first  (each  point  Qt  has  p  components).  It  is  done  using  the  Latin  square  or  Latin  hypercube  sampling,20 
which  is  useful  when  you  must  sample  a  ^-dimensional  space  exceedingly  sparsely,  at  M  points.  The  approach  is  to 
partition  each  (normalized)  design  parameter  (dimension)  into  M  segments,  so  that  the  whole  space  is  partitioned 
into  Af  cells.  The  M  cells  to  contain  the  sample  points  are  chosen  by  the  following  algorithm:  /)  randomly  choose 
one  of  the  Af  cells  for  the  first  point,  ii)  eliminate  all  cells  that  agree  with  this  point  on  any  of  its  parameters  (that  is, 
cross  out  all  cells  in  the  same  row,  column,  etc.),  leaving  (M-Yf  candidates,  iu)  randomly  choose  one  of  these 
remaining  candidates,  eliminate  new  rows  and  columns,  and  continue  the  process  until  there  is  only  one  cell  left, 
which  then  contains  the  final  sample  point.  The  result  of  this  construction  is  that  each  design  parameter  will  have 
been  tested  in  every  one  of  its  subranges.  Figure  3  provides  with  an  example  of  “wise”  population  of  two- 
dimensional  search  region  II  as  compared  to  that  of  “straightforward”  population  one  might  think  of. 

To  generate  the  original  sequence  of  points  q..,  i  =  1,...,M  ,  j  =  (where  qtj  is  the y-th  component  of  the 

z'-th  point  Qi)  the  PSI-method  employs  the  so-called  LPX  sequence  generation  procedure,21  which  in  turn  inherits 
Sobol’s  quasi-random  sequences22  generator  by  Antonov  and  Saleev.23  The  Sobol’  sequence  generates  quasirandom 
numbers  qtj  ?  i  =  1  ,...,M  ,  j  =  1, between  zero  and  one  directly  as  binary  fractions  of  length  w  bits,  from  a  set  of 

w  special  binary  fractions,  v*,  k  =  1,...,  w ,  called  direction  numbers.  In  Sobol’s  original  method,  the  z-th  number  Qf  is 
generated  by  XORing  (bitwise  exclusive  or)  together  the  set  of  v^’s  satisfying  the  criterion  on  k ,  “the  k- th  bit  of  i  is 
nonzero.”  In  other  words,  as  i  increments,  different  ones  of  the  vf  s  flash  in  and  out  of  q..  on  different  time  scales. 

By  construction,  the  first  direction  number  v\  alternates  between  being  present  and  absent  most  quickly,  while  vk 
goes  from  present  to  absent  (or  vice  versa)  only  every  2k~l  steps. 

The  advantage  of  Sobol’s  approach  (LPX  sequence  generation  procedure)  is  that  the  sequence  is  generated 
number-theoretically,  rather  than  randomly  (as  for  other  known  approaches24),  so  successive  points  at  any  stage  sort 
of  “know”  how  to  fill  in  the  gaps  in  the  previously  generated  distribution  and  keep  filling  them  in,  hierarchically 
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(see  example  on  Fig.4).  Figure  5  shows  that  the  coverage  of  LPX  sequence  is  also  much  more  equally  distributed  in 
comparison  to  the  embedded  Microsoft  Windows  random  numbers  generator  (RNG).) 

Finally,  the  unit  /7-dimensional  cube  is  stretched  to  the  parametric  constraints  (5)  by  following  scaling  procedure 

a)  =  a'j  +qij(a*  -a j) ,  j  =  1  i  =  .  (10) 


a2 


1 


ai 


Figure  3.  Example  of  straightforward  (a)  versus  “wise”  Latin-cube  (b)  sampling. 


Values  of  functional  dependencies  are  being  computed  for  these  M  trial  points.  If  they  satisfy  corresponding 
constraints  (6),  the  quality  criteria  ®z(cf),  /  =  l,...,v  are  also  being  calculated  at  each  trial  point  i  =  l,...,N, 
N  <M  . 


ii  M  Mi*  ^  MM  M»  v  e  M  u  J  >. 

a)  b) 

Figure  5.  Comparison  of  LPT  sequence  coverage  (a)  with  Windows  RNG  coverage  (b)  in  the  plain  of  two  (1st  vs.  10th)  out  of 

25  parameters  for  2048  trials. 
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The  parameter  space  is  investigated  in  three  stages.  First,  a  table  of  trials  is  ascribed  to  each  /- th  criterion 
Oz(a) ,  and  the  values  of  O^a1) ,...,  O^a^)  are  arranged  in  ascending  order  (assuming  that  all  the  criteria  must  be 
minimized). 

At  the  second  stage,  the  expert  chooses  preliminarily  the  criterial  constraints  ®**  (8).  During  the  third  stage  the 
problem’s  solvability  is  checked  meaning  that  the  set  of  all  a1  satisfying  all  inequalities  (9)  simultaneously  is 
determined.  If  the  set  of  these  vectors  a1  is  nonempty,  then  the  problem  of  the  feasible  set  construction  is  solvable. 
Otherwise,  one  has  to  either  correct  the  values  of  ®**  or  to  return  to  the  first  stage  and  increase  the  number  of  trials 
to  repeat  the  second  stage  with  a  larger  table.  The  procedure  is  continued  until  D  proves  to  be  nonempty  and  the 
maximum  values  of  ®**  are  specified.  After  that,  the  Edgeworth-Pareto  set  EPS  is  constructed  and  analyzed. 

V.  From  Multicriteria  Optimization  to  Multicriteria  Identification 

Obviously,  the  problem  formulation  in  Section  3  can  be  easily  adapted  to  the  multictriteria  parametrical 
identification  problems.  To  start  with  we  note  that  in  the  problem  of  multictriteria  parametrical  identification  or 
matching  experimental  data  to  the  predefined  mathematical  model  vector  of  variable  parameters  a  to  be  optimized 

may  include  t0  and  tj- . 

We  denote  by  CD"’ (a),  /  =  l,...,v  the  indices  (criteria)  resulting  from  the  analysis  of  the  mathematical  model 
that  can  be  represented  by  the  Eq.(4).  The  model  (4)  can  include  some  random  perturbations  like  white  noise  or  any 
other  disturbances  (inaccurate  wind  in  our  case).  In  this  case  we  will  consider  ®™(a)  being  a  mathematical 
expectation  of  corresponding  index. 

On  the  contrary,  let  ®zexp ,  /  =  l,...,v  denote  experimental  values  of  the  /-th  criterion  measured  on  the  prototype. 

Of  course  we  assume  the  experiment  to  be  sufficiently  accurate  and  complete  as  well  as  amount  of  measured  data 
available  to  be  sufficient  for  correct  formulation  on  the  identification  problem  at  hand.  If  the  data  for  several 
experiments  is  available  then  ®zexp  will  represent  the  mathematical  expectation  or  some  other  estimate  of  the  random 
variable. 

Now  instead  of  quality  or  performance  criteria  (7)  we  will  use  the  following  adequacy  (proximity,  closeness) 
criteria 

9?,(®r(a),®rp),  i=i,...,v,  (ii) 

where  51,  (®,"'(a),®,exp)  denotes  some  operator  applied  to  simulated  and  experimental  indices  (it  might  be  their 

ratio,  module  of  the  difference,  etc.). 

Therefore,  criterial  constraints  (9)  can  now  be  rewritten  as 

*,(®r(«),®rpE®r,  /= 02) 

For  this  problem  to  a  considerable  extent  the  values  of  ®**  depend  on  the  accuracy  of  the  experiment  and  physical 
sense  of  the  proximity  criteria  (11). 

This  brings  us  to  the  following  formulation  of  multicriteria  parameter  identification  problem  for  system  (4)  and 
sets  of  constraints  (5),  (6),  and  (12).  The  multicriteria  parameter  identification  problem  is  to  find  an  Edgeworth- 
Pareto  set  EPS,  EPS  cD,so  that: 

®/(EPS)  =  min9I/(®;z(a),®/exp),  /  =  l,...,v.  (13) 

aeD  '  > 

After  finding  EPS,  the  most  preferable  or  optimal  vector  a0 ,  a0  e  EPS  matching  the  physical  sense  of  the 

object  and/or  results  of  the  experiments  can  be  finally  determined  (chosen).  If  not,  then  the  problem  of  identification 
has  an  ambiguous  solution  (one  should  keep  in  mind  that  as  a  rule,  in  practice  some  of  the  criteria  are  calculated 
with  comparatively  high  accuracy,  while  others  are  determined  with  considerable  errors).  Theoretically,  to  resolve 
this  ambiguity  researcher  can  reconsider  the  rigidity  of  or  maybe  add  more  constraints.  Additional  experiments 
might  help  also.  However,  usually  this  can  be  done  on  rare  occasions  only,  because  basically  the  ambiguity  of 
restored  parameter  is  the  price  to  be  paid  for  incomplete  simulation  of  a  real  object  by  a  mathematical  model, 
incompleteness  of  the  full-scale  experiment,  etc. 
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VI.  Software  Description 

The  following  briefly  discusses  possibilities  of  the  usage  of  the  PSI  method  within  the  frame  of  the  MO VI  1.3 
software  package  adapted  to  Mathwork’s  Matlab/Simulink.  This  package  allows  user  to  perform  feasibility 
analysis/design  in  a  fairly  friendly  form.  The  total  number  of  variable  parameters  (they  maybe  both  continuous  and 
discrete)  exceeds  several  hundreds.  No  other  constraints  are  imposed  on  the  program.  Criteria  can  be  either 
minimized  or  maximized.  Some  or  even  all  of  the  criterial  constraints  if  unknown  a  priori  can  be  considered  as 
pseudocriteria. 

Figures  6-10  show  the  test  runs  of  the  identification  problem  where  several  rather  then  all  variable  parameters 
and  four  different  adequacy  criteria  were  used.  Fig. 6  demonstrates  an  example  of  test  tables  obtained  after  multiple 
runs  of  the  model  with  different  parameter  vectors.  Number  N  in  the  left-top  comer  indicates  the  total  number  of 
trials,  while  ND  (ND<N)  -  the  number  of  design  variable  vectors  in  the  feasible  set.  All  functional  failures  (trials 
that  did  not  meet  the  functional  constraints)  can  be  considered  separately  in  another  table.  By  softening  constraints, 
part  of  them  may  be  immediately  returned  to  the  feasible  set. 


Figure  6.  Example  of  fully  ordered  test  table. 


Histograms:  Feasible  Set 


Design  variable  |3-m1  j*]  {Show  Pareto-optimal  designs 


Ql 


950  960  970  980  990  1,000  1,010  1,020  1,030  1,040  1,050 

Design  Variable  3  (9.500000000E+02  -  1  050000000E+03) 


Figure  7.  Feasible  set  histogram  for  the  single  parameter. 


The  minimum  and  maximum  numbers  of  each  adequacy  criterion  are  presented  at  the  title  of  each  table  (vertical 
column).  On  this  step  MOVI  software  allows  an  expert  to  tmncate  the  whole  table  achieving  better  results  by 
working  with  the  small  portion  of  it  at  a  time,  and  to  correct  the  value  of  any  criterial  constraint  to  narrow/broaden 
the  feasible  set  (it  can  be  done  for  all  columns  from  the  left  to  the  right  decreasing  ND  value  for  every  criterion). 
Finally,  resulting  table  may  be  converted  to  the  simplified  form  containing  information  on  the  subsets  of  feasible 
solutions  and  Pareto-optimal  solutions  (NP<ND). 

Further  analysis  involves  graphical  representation  of  the  data  (see  examples  on  Figs. 7-10).  Fig. 7  represents  a 
histogram  for  the  specific  parameter  (how  many  of  the  trials  fall  into  the  certain  range).  Fig. 8  depicts  an  example  of 
criteria  versus  criteria  graph.  Figure  9  shows  the  criteria  versus  single  parameter  plot  for  all  trial  points  (meaning 
that  each  point  corresponds  to  a  single  parameter  vector).  In  addition  the  criteria  versus  single  parameter  plot  can  be 
obtained  for  any  specified  parameter  by  running  some  additional  runs  with  all  other  parameters  fixed  (as  it  shown  on 
Fig.  10).  It  is  possible  to  change  ranges  for  the  chosen  parameter  here.  Moreover,  the  values  of  any  other  component 
of  the  parameter  vector  can  be  corrected  also. 


Figure  8.  Criterion  vs.  criterion  graph. 


Figure  9.  Criterion  vs.  parameter  graph  constructed  using 
all  feasible  solutions. 
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Figure  10.  Criterion  vs.  parameter  graph. 

VII.  Discussion  of  Parametric  ID  Results 

As  mentioned  in  Section  2,  over  30  parameters  were  used  along  with  eight  different  criteria.  Among  them  there 
are  two  criteria  describing  the  closeness  of  the  horizontal  and  vertical  projections  of  the  trajectories,  and  three 
criteria  relating  to  the  adequacy  of  the  natural  eigenvalues  (power  spectrum)  for  all  channels  (roll,  pitch  and  yaw). 

Figure  11  shows  the  real  drop  trajectory,  prototype  trajectory  (all  gain  factors  are  equal  to  unity)  and  the 
trajectories  that  were  found  to  be  the  best  with  respect  to  the  each  specific  criteria  named  above. 


Figure  11.  Optimal  parafoil  trajectories  with  respect  to  different  criteria. 

The  choice  was  made  among  the  results  of  several  thousand  trials  (about  two  months  of  continuous  PC  run). 
However  the  set  of  feasible  (and  Pareto-optimal)  solutions  included  much  fewer  trials  because  of  criterial  constraints 
applied  to  some  of  the  state  variables  (angle  of  attack,  speed  components,  altitude).  As  seen  from  Fig.  11,  all 
trajectories  are  fairly  close  to  each  other  which  actually  attests  the  high  quality  of  the  model.  As  expected  and 
predicted  by  others,25  not  much  difference  was  observed  compared  to  the  6-DoF  model  the  authors  developed 
earlier.5  The  only  noticeable  difference  was  in  slightly  smaller  discrepancy  in  the  natural  eigenvalues  (power 
spectrum)  for  pitch  and  yaw  channels.  Of  course  the  8  DoF  model  exhibited  closer  match  to  that  of  real  drop  data. 
Analysis  of  the  power  spectrum  exposed  which  frequencies  were  missing  and  therefore  defined  limitation  of  the  6- 
DoF  versus  8-DoF  model.  It  also  allowed  evaluating  the  accuracy  and  applicability  of  the  wind  data  gathered  by  the 
dropsonde. 

Yet,  the  trajectories  on  Fig.  11  do  not  match  completely.  Of  course,  neither  the  values  of  variable  parameters 
match  well.  Moreover,  even  if  only  cost  function  is  considered  (and  maybe  only  a  single-criterion  ID  method  is 
applied),  multiple  near-optimal  (local)  variable  data  sets  can  be  found  as  seen  from  Fig.  12.  This  figure  presents  the 
values  for  33  variable  parameters  (coefficients  gains)  for  several  sets  for  the  certain  cost  function  the  value  of  which 
for  each  set  is  also  shown  on  the  bottom.  Fig.  13  represents  the  same  data  but  graphically  to  show  the  variation  of 
parameters  with  respect  to  their  nominal  unity  values  (after  several  preliminary  runs  the  range  for  all  gain  factors 
was  established  as  [0.2;2]). 

While  several  preliminary  trials  when  only  a  few  parameters  were  allowed  to  be  changed  (corresponding  to  set 
1-5  columns  in  the  table  of  Fig.  12)  decreased  the  cost  function  from  10.44  dimensionless  units  to  about  4  units, 
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further  adjustment  involving  over  thirty  variable  parameters  was  needed  to  decrease  the  value  of  the  cost  function 
further.  Some  of  the  resulting  local  optimum  sets  are  shown  in  the  last  five  columns  of  Fig.  12.  What’s  interesting  is 
that  having  approximately  the  same  value  of  the  cost  function  (less  than  2  units)  and  resulting  in  fairly  close 
trajectories  these  sets  differ  from  each  other  (sometime  as  much  as  twofold).  That  means  that  having  so  many 
variable  parameters  we  can  redistribute  them  in  several  ways  to  achieve  approximately  the  same  magnitude  of  the 
cost  function. 


Parameter/Value 

SetO 

Set  1 

Set  2 

Set  3 

Set  4 

Set  5 

Set  6 

Set  6a 

Set  7a 

Set  33-1 

Set  33-2 

cl 

deicon 

1 

1 

1 

0.23983 

0.23983 

0.23983 

0.23983 

0.2878 

0.2865 

0.2865 

0.494613 

0.349706 

c2 

xcq 

0 

0 

0 

0 

0 

0.386 

0.386 

0.4632 

0.59664 

0.59664 

0.546635 

0.535231 

c3 

AlphaO 

10 

1 

1 

1 

1 

1 

0.936 

1.1392 

1.1739 

1.1739 

1.240674 

1.255886 

c4 

CLalpha 

0.0375 

1 

1 

1 

1 

1 

1.0022 

0.91294 

0.96735 

0.96735 

0.927995 

0.941119 

c5 

CLDf 

0.2 

1 

1 

1 

1 

1 

1.07 

1.137 

1.2707 

1.2707 

1.770711 

1.769055 

c6 

CD0 

0.14 

1 

1 

1 

1 

1 

1.062 

0.86667 

1.0494 

1.0494 

0.96331 

0.961487 

c7 

A2 

0.25 

1 

1 

1 

1 

1 

0.995 

0.80066 

0.4483 

0.4483 

0.634647 

0.629522 

c8 

ADf 

0.2 

1 

1 

1 

1 

1 

1.064 

1.2 

0.75852 

0.75852 

0.413619 

0.406355 

c9 

CmO 

-0.33 

1 

1 

1 

1 

1 

1.0486 

0.8 

1.1106 

1.1106 

0.937026 

0.940533 

clO 

Rig 

8 

1 

1 

1 

1 

1 

1 

1.0523 

0.58525 

0.58525 

0.958859 

0.957133 

ell 

Cmq 

-6.39 

1 

1 

1 

1 

1 

0.894 

0.8 

1.2255 

1.2255 

0.588472 

0.583637 

cl  2 

CYbnom 

-0.005 

1 

1 

1 

1 

1 

1 

0.8 

1.1803 

1.1803 

0.461411 

0.442807 

cl  3 

CYbqrad 

-0.0001 

1 

1 

1 

1 

1 

1 

1.2 

0.86021 

0.86021 

1.440569 

1.451131 

cl  4 

CYDanom 

-0.007 

1 

1 

1 

1 

1 

1 

1.2 

0.2263 

0.2263 

1.132878 

1.282393 

cl  5 

CYDagrad 

0.0012 

1 

1 

1 

1 

1 

1 

1.2 

1.4006 

1.4006 

1.214801 

1.340795 

cl  6 

Clbnom 

-0.0014 

1 

1 

1 

1 

1 

1 

0.95664 

1.3619 

1.3668 

0.456535 

0.462645 

cl  7 

Clbqrad 

-0.001 

0 

0 

0 

0 

0 

0 

0 

0 

0.000161 

0.00013 

8.98E-05 

cl  8 

CIDanom 

-0.0063 

1 

1 

1 

1 

1 

1 

1.2 

0.74132 

0.74518 

0.530218 

0.62248 

cl  9 

CIDagrad 

-0.001 

0 

0 

0 

0 

0 

0 

0 

0 

0.00011 

-4.51  E-06 

-3.78E-06 

c20 

Clp 

-0.15 

1 

1 

1 

1.0765 

1.0765 

1.0765 

0.8612 

1.3461 

1.3476 

1.060879 

1.07091 

c21 

Clr 

0.0775 

1 

1 

1 

1 

1 

1.0037 

0.80296 

0.68709 

0.68807 

0.56154 

0.565489 

c22 

Cnbnom 

0.007 

1 

1 

1 

1 

1 

1 

0.8 

0.57465 

0.57465 

0.473491 

0.469812 

c23 

Cnbqrad 

-0.0003 

1 

1 

1 

1 

1 

1 

1.2 

1.0893 

1.0893 

0.77051 

0.774133 

c24 

CnDanom 

0.03 

1 

0.94344 

0.94344 

0.94344 

0.94344 

0.94344 

0.98673 

0.89177 

0.89177 

0.963727 

0.965694 

c25 

CnDaqrad 

-0.001 

1 

0.50087 

0.50087 

0.50087 

0.50087 

0.50087 

0.4007 

0.34391 

0.34391 

0.644465 

0.540929 

c26 

Cnp 

0.023 

1 

1 

1 

1 

1 

0.99 

1.188 

1.519 

1.519 

1.024963 

1.022094 

c27 

Cnr 

-0.0936 

1 

1 

1 

1 

1 

1.0012 

1.2014 

1.155 

1.155 

1.10557 

1.108614 

c28 

kA 

0.899 

1 

1 

1 

1 

1 

1 

1 

1 

1 

0.539401 

0.548065 

c29 

kB 

0.339 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1.382761 

1.384552 

c30 

kC 

0.783 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1.259202 

1.278601 

c31 

kA* 

0.63 

1.43 

1.43 

1.43 

1.43 

1.43 

1.43 

1.43 

1.43 

1.43 

1.798577 

1.760393 

c32 

kB* 

0.817 

0.41 

0.41 

0.41 

0.41 

0.41 

0.41 

0.41 

0.41 

0.41 

0.497986 

0.525471 

c33 

kC* 

1 

0.78 

0.78 

0.78 

0.78 

0.78 

0.78 

0.78 

0.78 

0.78 

0.402627 

0.202126 

Cost  Function 

10,444 

4,224 

4,058 

3,947 

3,883 

3,700 

1,812 

1,271 

1,271 

1,272 

1,181 

Figure  12.  “Quasi-optimal”  values  for  33  variable  parameters  for  different  sets  of  variable  parameters  for  a  single 

criteria. 


1,272  1,271  a  1,271  o  1,181  +1,812  o3,700  *3,883  3,947  ±4,058  *  4,224  <>10,444 

- U - ' - 1 - 1 - 1 - 


1.5 


A 


O 

$ 


+  +  + 


+  1 
o 


^  a  °  o  +  !  ® 

o - ogDgDDOooaocsaD - o - 

°  °  “  !  A 


+  A 

+  A  +  !  +  + 


0.5 


^  O 


+  + 

A 

O  A 
O 


O 

a 


a  — a 


5  10  15  20  25 

Coefficient 

Figure  13.  Graphical  representation  of  parameters  spread. 
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On  the  one  hand,  detailed  analysis  of  solutions  revealed  explicit  and  implicit  correlation  between  some  variable 
parameters  (darkened  cells  on  Fig.  14  show  correlation  between  some  of  25  variable  parameters).  That  allows 
decreasing  the  original  number  of  variable  parameters  by  about  30%  and  alternatively  (if  needed)  adding  new 
parameters  that  were  not  considered  in  the  original  setup  without  increasing  dimension  of  the  problem. 
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Figure  14.  Strong  correlation  between  some  of  variable  parameters. 


This  also  provides  insight  into  the  set  of  additional  states  of  the  system  that  need  to  be  observed  and  recorded. 
The  data  on  system’s  motion  is  usually  obtained  either  with  the  help  of  GPS/IMU  unit  atop  the  payload  or  by 
tracking  the  payload  using  a  cinetheodolite  system  with  the  following  post  analysis.26  As  shown  by  our  previous 
study4,5  having  this  data  is  almost  enough  for  validation  of  6-DoF  models.  For  a  more  fundamental  study  including 
separate  payload  and  parafoil  behavior  (higher  fidelity  models)  and  perhaps  inflation  dynamics,  the  canopy  motion 
has  to  be  investigated  separately.  Otherwise  due  to  lack  of  experimental  studies  and  measurements,  the 
parafoil/payload  interactions  are  often  postulated  in  analytical  modeling,  resulting  in  theoretical  predictions  based  on 
uncertain  assumptions.  Therefore,  there  is  a  need  to  experimentally  investigate  motion  of  the  parafoil,  employing  for 
instance  the  technique  of  measuring  two  angles  defining  a  direction  of  the  riser  with  respect  to  the  payload,  offered 
by  Lee  et  al.27,  video  imaging  of  the  canopy  from  the  camera  installed  atop  payload,28,29  or  by  applying  the 
algorithms  originally  developed  to  for  payload26  to  estimate  the  pose  of  the  parafoil. 

VIII.  Conclusion 

Discussion  presented  in  this  paper  persecutes  the  goal  of  correctly  formulating  the  problem  of  multicretiria 
parametrical  identification  of  the  parafoil-based  delivery  system.  It  is  suggested  to  implement  the  well-established 
PSI  multicriteria  optimization  method  to  investigate  the  set  of  feasible  parameters  and  solve  the  identification 
problem.  The  paper  shows  that  in  total  the  problem  contains  as  much  as  several  dozens  variable  parameters  and 
several  distinctive  adequacy  criteria.  Different  sets  of  parameters  affect  these  criteria  non-adequately.  Moreover, 
minimum-criterion  solutions  do  not  coincide.  Currently  authors  are  performing  more  simulations  with  existing  set  of 
flight  test  data  and  expect  to  involve  some  more  to  be  able  to  complete  tweaking  the  model. 
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